package geos

import (
	"github.com/paulmach/orb/encoding/wkt"
	"github.com/paulmach/orb/planar"
	"github.com/paulmach/orb/project"
	"math"
)

func rad(d float64) (r float64) {
	r = d * math.Pi / 180.0
	return
}

const EARTH_RADIUS = 6378137

const POLAR_RADIUS = 6356725

func Distance(lon1, lat1, lon2, lat2 float64) (distance float64) {
	//赤道半径(单位m)
	rad_lat1 := rad(lat1)
	rad_lon1 := rad(lon1)
	rad_lat2 := rad(lat2)
	rad_lon2 := rad(lon2)
	if rad_lat1 < 0 {
		rad_lat1 = math.Pi/2 + math.Abs(rad_lat1)
	}
	if rad_lat1 > 0 {
		rad_lat1 = math.Pi/2 - math.Abs(rad_lat1)
	}
	if rad_lon1 < 0 {
		rad_lon1 = math.Pi*2 - math.Abs(rad_lon1)
	}
	if rad_lat2 < 0 {
		rad_lat2 = math.Pi/2 + math.Abs(rad_lat2)
	}
	if rad_lat2 > 0 {
		rad_lat2 = math.Pi/2 - math.Abs(rad_lat2)
	}
	if rad_lon2 < 0 {
		rad_lon2 = math.Pi*2 - math.Abs(rad_lon2)
	}
	x1 := EARTH_RADIUS * math.Cos(rad_lon1) * math.Sin(rad_lat1)
	y1 := EARTH_RADIUS * math.Sin(rad_lon1) * math.Sin(rad_lat1)
	z1 := EARTH_RADIUS * math.Cos(rad_lat1)

	x2 := EARTH_RADIUS * math.Cos(rad_lon2) * math.Sin(rad_lat2)
	y2 := EARTH_RADIUS * math.Sin(rad_lon2) * math.Sin(rad_lat2)
	z2 := EARTH_RADIUS * math.Cos(rad_lat2)
	d := math.Sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2))
	theta := math.Acos((EARTH_RADIUS*EARTH_RADIUS + EARTH_RADIUS*EARTH_RADIUS - d*d) / (2 * EARTH_RADIUS * EARTH_RADIUS))
	distance = theta * EARTH_RADIUS
	return
}

/***
 *  这个根据一个经纬度坐标、距离然后求另外一个经纬度坐标的函数
 *  例如，我要找一个坐标点(lat,lon)的5公里范围内的所有商户信息、景点信息等。这个MBR就是一个最大的范围，
 *  这个矩形是包含5公里范围内所有这些有效信息的一个最小矩形。
 *  求出四个方向0度、90度、180度、270度方向上的四个坐标点就可以得到这个矩形。
 * @param lat
 * @param lon
 * @param distance 距离，单位：公里
 * @param angle 正北极方向顺时针的角度
 * @return
 */
func AddDistance(lat, lon, distance, angle float64) (newLat, newLon float64) {
	dx := distance * 1000 * math.Sin(angle*math.Pi/180.0)
	dy := distance * 1000 * math.Cos(angle*math.Pi/180.0)
	ec := POLAR_RADIUS + (EARTH_RADIUS-POLAR_RADIUS)*(90.0-lat)/90.0
	ed := ec * math.Cos(lat*math.Pi/180)
	newLon = (dx/ed + lon*math.Pi/180.0) * 180.0 / math.Pi
	newLat = (dy/ec + lat*math.Pi/180.0) * 180.0 / math.Pi
	return
}

// 获取点到线的距离(单位米)
func PointToLineStringDistance(lineStringStr string, pointStr string) (float64, error) {
	lineString, err := wkt.UnmarshalLineString(lineStringStr)
	if err != nil {
		return 0.0, err
	}
	point, err := wkt.UnmarshalPoint(pointStr)
	if err != nil {
		return 0.0, err
	}
	lineStringM := project.LineString(lineString, project.WGS84.ToMercator)
	pointM := project.Point(point, project.WGS84.ToMercator)
	distanceFrom := planar.DistanceFrom(lineStringM, pointM)
	return distanceFrom, nil
}

// 获取点到多边形的距离(单位米)
func PointToPolygonDistance(polygonStr string, pointStr string) (float64, error) {
	polygon, err := wkt.UnmarshalPolygon(polygonStr)
	if err != nil {
		return 0.0, err
	}
	point, err := wkt.UnmarshalPoint(pointStr)
	if err != nil {
		return 0.0, err
	}
	polygonM := project.Polygon(polygon, project.WGS84.ToMercator)
	pointM := project.Point(point, project.WGS84.ToMercator)
	distanceFrom := planar.DistanceFrom(polygonM, pointM)
	return distanceFrom, nil
}
